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1.  Introduction .  A  fundamental  problem  in  the  numerical  solution 


of  hyperbolic  equations  is  the  proper  approximation  of  the  boundary  condi¬ 
tions.  For  example,  the  leapfrog  scheme  applied  to  the  following  hyper¬ 
bolic  problem  is  unstable. 

WO  0<X<7T 

u(0,t)=sin(-t)  0<t  (1) 

u(x,0)=sin (x) 

The  mesh  for  this  scheme  is  x.=jAx=j7r/J  for  0<j<J.  The  exact  solution 

J 

is  sin  (x-t).  If  a  second  order  difference  approximation  for  the  spatial 
derivative  is  combined  with  a  leapfrog  scheme  for  time,  then  the  following 
scheme  is  obtained 

Ug=sin(-tn) 

,,n+T  n-1  9a.  /  n  n  x  /A 

UJ  ~uj  t( u j“U j_-j  )! Ax 

,  n+l_  n-1  .  f  n  n  ,  .  n  .  , 

Uj  -Uj  -At(Uj+-|-U  ._i  )/AX 


This  scheme  is  unstable.  If  the  outflow  boundary  is  modified  as  indicated 
below,  then  the  scheme  is  stable  and  has  second  order  accuracy. 

u3+1“uj-4t(“3-u3-i>/Ax 


When  the  method  of  lines  is  used  for  the  simple  linear  hyperbolic 
equation  (1)  with  periodic  boundary  equations,  then  the  resulting  differ¬ 
ence  scheme  is  stable,  provided  an  ODE  solver  with  automatic  step-size 
adjustment  is  used  to  solve  the  system  of  ordinary  differential  equations. 
Even  if  the  ODE  solver  uses  an  Euler  forward  time-step  scheme,  the  inte¬ 
gration  will  converge  as  the  mesh  spacing  is  taken  to  zero,  since  the 
ODE  solver  will  take  the  step  time-step  small  enough  as  a  function  of 
the  mesh  size  to  guarantee  convergence.  It  will  not  be  the  case  that 
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At=0 ( Ax ) .  Note  that  the  semi -discrete  approximation  produced  by  the 
method  of  lines  with  periodic  boundary  conditions  can  be  written  in 
the  form 

U.=AjJ  U_= (u_js  •  •  •  >Uj_-j  )  X  .=j7r/J 
where  the  matrix  A  is 

0  10... 

-1  0  10.. 

0-1  01.. 

1 

The  solution  of  this  differential  equation  is  given  in  terms  of  an 
exponential  matrix  as 

u(t)=u(0)eAt 

Since  the  matrix  is  skew  symmetric  and  cyclic  its  eigenvalues  are  pure 
imaginary  and  its  eigenvectors  are  orthogonal.  Therefore,  the  solution 
is  bounded  independently  of  the  number  of  mesh  points.  This  implies 
that  the  solution  of  this  semi -discrete  approximation  will  converge  to 
the  solution  of  the  original  equation  (1).  Therefore  any  spatial 
discretization  which  yields  a  skew  symmetric,  cyclic  matrix  will  define 
a  convergent  method  of  lines  approximation.  Stability  in  a 
finite  difference  scheme  for  hyperbolic  problems  is  in  a  sense  associated 
with  the  temporal  discretization. 

Unfortunately,  the  method  of  lines  does  not  necessarily  produce 
a  stable  scheme  when  the  boundary  conditions  are  not  periodic.  However, 
the  method  of  lines  does  seem  to  be  more  likely  to  yield  a  stable  scheme 
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than  a  leapfrog  time  discretization. 

Our  purpose  is  an  experimental  study  of  some  boundary  difference 
approximations  for  use  on  hyperbolic  systems  where  the  method  of  lines 
is  used  for  the  temporal  discretization.  Our  results  will  refer  mainly 
to  the  Runge-Kutta-Fehl berg  ODE  solver,  although  we  intend  to  experi¬ 
ment  with  the  Adams  method  of  Shampine  [9]  in  the  future.  We  have  found 
it  important  to  include  test  cases  for  hyperbolic  systems  (more  than 
one  independent  variable)  for  which  the  characteristics  lie  on  both  sides 
of  the  boundary.  This  is  in  agreement  with  comments  by  Chu  [3]  and 
Sundstrom  [10].  We  are  interested  in  boundary  approximations  which  can 
be  incorporated  into  a  general  PDE  solver  to  treat  hyperbolic  systems 
in  two  dimensions.  Such  solvers  for  parabolic  equations  in  one  dimen¬ 
sion  are  described  by  Sincovec  and  Madsen  [11],  Carver  [2],  Loeb  [6], 

Bowen  [1],  Hastings  [12]  and  others.  Because  of  our  interest  in 
general  hyperbolic  systems,  we  cannot  consider  boundary  approximations 
stated  in  terms  of  specific  variables  for  specific  equations.  We  can 
only  consider  algorithms  which  can  be  presented  in  a  general  gramework. 

We  will  test  two  such  algorithms. 

Of  course,  such  a  general  algorithm  requires  the  user  to  apply 
it  in  such  a  way  as  to  produce  a  properly  posed  hyperbolic  problem. 

We  must  allow  the  user  the  flexibility  to  set  the  boundary  conditions. 
Eventually,  we  might  be  able  to  supply  an  optional  check  to  see  that 
the  boundary  conditions  are  consistent  with  the  hyperbolic  system. 

2.  Computational  results  indicating  stability  and  accuracy  of 
the  method  of  lines.  In  this  section  we  consider  difference  approximations 
for  the  system  (1).  These  are  semi -di screte  approximations  of  the  form 
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u_'=Au_+£  (2) 

T 

where  u_  =  u_(t)  =  ( ..  .u  .(t) ... )  is  a  vector  of  mesh  point  val  ues.  In 

vj 

this  section  we  will  look  at  the  eigenvalues  of  A  and  the  norm  of  the 
exponential  matrix 


for  four  finite  difference  approximations.  If  this  n*orm  is  bounded 
independent  of  the  spatial  mesh,  then  the  semi -discrete  approximation 
is  stable.  This  follows  from  the  integral  form  of  the  solution  of  (2) 


u(t)  =  u(o)e^t+ 


•0 


(4) 


If  the  eigenvectors  of  A  are  orthogonal,  then  a  bound  for  the  norm  of 
the  exponential  matrix  can  be  obtained  from  the  eigenvalues  of  A. 
Therefore,  we  compute  these  eigenvalues  and  also  the  norm  (3)  in  order 
to  gain  insight  into  the  stability  of  the  following  four  schemes. 


A.  An  inconsistent  scheme.  Here  a  one-sided  difference  is  used 
at  both  boundaries  in  spite  of  the  fact  that  the  solution  should  be 
specified  at  the  inflow  or  left  boundary.  This  must  yield  an  unstable 
approximation.  The  approximation  is  consistent,  and  if  it  were  also 
stable,  then  it  would  be  convergent.  That  is,  if  the  norm  of  the 
exponential  matrix 
eAt 

were  bounded  independently  of  the  mesh  spacing  the  approximation  would 
be  convergent,  which  is  impossible  since  no  boundary  condition  has  been 
specified  on  the  inflow  boundary.  The  scheme  is 
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(5) 

!<j<J 

B.  A  second  order  scheme.  This  scheme  is  the  same  as  the  previous 
one  except 

uQ(t)=sin(-t) 

It  is  only  first  order  at  the  boundary,  but  the  overall  accuracy  should 
be  second  order. 

C.  Fourth  order  with  a  third  order  boundary.  This  scheme  is 
given  below,  Oliger  £7]  has  shown  that  subtle  changes  are  required  in 
this  spatial  approximation  when  it  is  used  with  a  leapfrog  time  dis¬ 
cretization,  in  order  that  the  resultant  difference  scheme  be  stable. 
However,  it  seems  to  be  stable  without  modification  when  it  is  used 
with  a  variable  step  ODE  solver. 

Ug(t)=sin(-t) 

u.j  (t)=-(2uQ-3u1+6u3-u4)/6Ax 

f  ( g ) 

uj(t)="(2uj-2'16uj-i+16Vr2V2)/{24Ax) 

uj_l  (t)=-(uj_3"6Uj_2+3Uj_.j+2Uj)/(6Ax) 

uj(t)=-{-2uJ_3+9uJ_2-18UJ_1+nuJ)/(6Ax) 


u0(t)=~(Ul(t)-Uo(t)) 

uj  (t)=:2S{uj+i  W'Uj-T 
uj(tK  4^uj(t)-u  (t» 
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D •  A  fourth  order  scheme  with  a  fourth  order  boundary  approximation. 
This  is  the  same  scheme  as  the  one  above  except  that  one-sided  fourth 
order  differences  are  used  at  the  boundary. 

,  u1(t)=-(-6uQ-20u]+36u2-l2U  +2u  )/(24ax) 
uj-l(t)=-(-2uj.4+'|2uJ.3-36uJ_2+20uJ_1+6UJ)/(24Ax)  (7) 

u^t)=-(6UJ_4-32yj_3+72uJ_2-96uJ_1+50uJ)/(24Ax) 

The  above  four  schemes  can  all  be  written  in  the  matrix  form  of 
equation  (2).  The  maximum  of  the  real  parts  of  the  eigenvalues  of 
the  matrix  A  for  these  four  schemes  are  given  in  Table  I.  The  eigenvalues 
of  A  were  determined  by  using  the  IMSL  QR  routine  on  the  CDC  6400 
at  the  University  of  Colorado.  The  exponential  matrix  was  determined 
by  summing  its  series  expansion.  The  norm  is  that  induced  by  the  vector 
maximum  norm.  For  the  inconsistent  scheme  (A)  the  eigenvalues  are  all 
pure  imaginary  with  a  double  or  triple  root  at  zero  depending  on  whether 
J  is  even  or  odd.  The  instability  of  this  scheme  is  evident  from  the 
norm  of  the  exponential  matrix  but  not  from  the  eigenvalues. 

Schemes  (B)  and  (C)  would  appear  to  be  stable  from  this  analysis. 
However,  we  might  expect  the  solution  of  scheme  (D)  to  show  exponential 
growth  in  time  since  its  matrix  has  an  eigenvalue  with  positive  real 
part. 

In  order  to  provide  a  more  complete  test  of  these  four  schemes 
we  wrote  a  code  for  these  schemes  applied  to  equation  (1).  This  pro¬ 
vides  a  direct  test  of  the  stability  and  accuracy  of  these  schemes. 

Table  II  shows  the  error  obtained  with  the  various  schemes  after 
integration  to  the  indicated  value  of  t=T  using  the  mesh  resolution 
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determined  by  J,  Note  that  the  number  of  intervals  per  wave  is  2 ( J-l ) 
since  the  mesh  runs  from  x=0  to  x=tt  and  J+l  is  the  number  of  mesh 
points.  Scheme  (A)  is  clearly  unstable.  Schemes  (B)  and  (C)  seem  to 
be  stable  which  is  consistent  with  the  results  in  Table  I  giving  the 
characteristics  of  the  matrices  corresponding  to  these  schemes.  Scheme 
(D)  seems  to  be  weakly  unstable  when  the  system  is  solved  with  the  RKF 
ODE  solver.  However  this  scheme  seems  to  be  stable  when  the  Runge-Kutta 
scheme  with  a  fixed  ratio  At/ Ax  is  used. 

3.  A  variable  coefficient  problem.  A  hyperbolic  problem  which 
is  more  typical  of  many  applications  than  equation  (1)  is  the  following 
defined  on  the  interval  0<x£tt. 

ut+cos(t)ux=cos(x-t)(cos(t)-l)=r(x,t) 

If  cos (t)_>0  then  u(Qst)=sin(-t)  (8) 

If  cos(t)<0  then  u(Trst)=sin  (ir-t) 
u(x,0)=sin  x 


The  solution  of  this  problem  is  u(x 9t)=sin ( x-t ) .  The  mesh  is  x  .= j tt/ J 9 
for  0<Jj<J .  In  this  problem  the  inflow  and  outflow  boundary  alternate 
between  the  two  endpoints  of  the  interval.  When  cos(t)>_0  the  left 
boundary  is  the  inflow  point.  This  makes  the  use  of  an  ODE  solver 
awkward  if  the  method  of  equation  (6)  is  used  to  define  the  system  of 
differential  equations.  When  cos(t)>0  the  unknowns  are  (u-j  (t) ,. . .  9Uj(t)) 
and  when  cos(t)<0  the  unknown  vector  has  shifted  to  (uQ(t) , . . .  ,Uj_-j  (t) ) . 
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Therefore  we  differentiate  the  boundary  condition  so  that  the  system 
of  differential  equations  always  contains  the  same  unknowns. 

E.  A  second  order  scheme  for  equation  (8). 

If  cos(t)>_0  then 

uo(t)=^t(sl*n("t)  )="cos(t) 
otherwise 

Uq ( t )  =  - ( u-j  -Uq)/ Ax+r(0,t)  (9) 

If  c o s ( t ) <0  then 

Uj(t)=^(sin(7r-t))=COS(t) 

otherwise 

Uj(t)  =  -(Uj-Uj_i  )/Ax+r(Trst) 

This  scheme  uses  a  differentiated  form  of  the  boundary  condition  at  an 
inflow  boundary  and  a  one-sided  first  order  difference  approximation 
to  the  differential  equation  at  an  outflow  boundary.  The  definition 
of  the  differential  equation  used  to  define  the  solution  along  the 
boundary  line  varies  depending  on  the  inflow-outflow  nature  of  the 
boundary.  However,  the  solution  along  these  boundary  lines  is  always 
determined  by  a  differential  equation. 

F.  A  fourth  order  scheme  for  equation  (8) . 

If  cos(t)>Q  then 

UQ(t)=-COS(t) 

otherwise 

ug(t)  =  -6'3(u)0+r(05t) 
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Here  S'  is  the  third  order  difference  approximation  of  u  (0)  using 
o  x 

(xQ,x1 ,x2,x3) .  The  equation  for  uj(t)  is  similar.  The  remainder  of 
the  system  is  identical  with  that  of  equation  (6). 

The  results  of  using  these  schemes  to  approximate  the  solution 
of  equation  (8)  is  given  in  Table  III.  These  results  indicate  that 
these  schemes  are  stable.  The  norm  of  the  second  order  approximation 
shows  a  slow  linear  growth  with  time.  The  error  shows  the  expected 
asymptotic  behavior  with  J  (approximately) .  The  behavior  of  this  method 
on  a  more  complex  multidimensional  problem  awaits  testing  which  we 
hope  to  carry  out  in  the  near  future. 


4.  General  boundary  approximation  algorithms.  In  this  section 
we  consider  a  program  for  the  following,  more  general  class  of  nonlinear 
hyperbolic  equations. 


ft  =  f^(a(u,x,t))+h(u,x,t)  (11) 

or  the  nonconservation  form 

=  f(|^5U_5X,t)+h(u_5X5t)  (11) 


Here  L>  and  h  are  general  vector  valued  functions  and  u_( x , t )  is  the 
vector  solution.  We  assume  that  boundary  conditions  are  given  at  two 
end  points  xs9  and  x=b.  We  consider  two  methods  to  specify  these  boundary 
conditions . 

The  first  method  requires  the  specification  of  a  subset  of  the 
unknowns  at  each  boundary  point.  Consider  the  left  boundary  x=a.  The 
unknowns  are  (uT(x,t) , . . . ,uM(x,t)) .  The  p  unknowns  (u  ,...,u  )  from 

*  1 1 '  m-j  nip 

the  set  I  =  (Pt . pm}  are  specified  as  follows: 

s  m 
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(  9  J  t  )  S^i  (  U  J  J  (  9  5  t  )  ,  t  j 


(a,t)  =  Sp(un 


(a,t) ,t) 


(13) 


Here  u_TT  =  (u  , . . .  su  )  is  the  compliment  of  uT  =  (um  ,...,um  ).  The 
11  ml  mM-p  1  ml  mp 

problem  specification  must  include  the  integer  p  and  the  functions 

Sr ...,S  at  both  boundary  points.  Note  that  p  may  depend  on  the  time  t. 

The  functions  S.  are  used  to  set  the  values  of  Uj  at  the  boundary.  The 

values  of  u ^ ^  are  computed  from  the  hyperbolic  equation,  using  one  sided 

approximations  for  spatial  derivatives. 

For  example,  consider  the  variable  coefficient  problem  given  by 

equation  (8).  At  the  left  boundary  (x=0),  if  cos(t)>0,  then  for  the 

number  of  boundary  constraints  we  have  p=l  .  The  function  (uj j ,t)  = 

S-j  (t)=-sin(t) .  Note  that  Ujj  is  empty  in  this  case.  If  cos(t)<0, 

then  p=0  at  the  left  boundary  and  Uj  is  empty.  In  this  case  the  value 

of  tL(t)  (here  U.  (t)  denotes  the  approximation  to  u(x.,t)  on  the  "time 

o  J  J 

line")  is  obtained  from  the  differential  equation 
dU0 

jjr-  =  -cos(t)53(U)0+r(x0,t)  (14) 

where  6^  represents  the  onesided  difference  approximation. 

When  the  ODE  solver,  such  as  the  Runge-Kutta-Fehl berg  is  used, 
there  is  a  slight  problem  in  implementing  this  algorithm.  When  the 
character!' Stic  slope  cos(t)  changes  sign,  the  nature  of  the  system  of 
ordinary  differential  equations  changes.  When  cos(t)>0,  the  unknowns  in 
the  system  (8)  are  (U-j , . . .  ,U.) ,  but  when  cos(t)<0  the  unknowns  are  (Ug, . . .  ) . 
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The  ODE  solver  always  works  with  the  full  set  of  unknowns  including  the 
boundary  values,  that  is  (Ug,...,Uj).  However,  in  computing  the  "right 
side"  functions  in  the  Runge-Kutta  steps  the  boundary  constraints  are 
applied  to  set  the  boundary  values  for  variables  in  the  Uj  sets.  If  the 
system  of  ODE's  is  written 


Ej(Uq, . . . ,Uj,t) 


05) 


and  Ug  is  in  the  constrained  set  Uj  for  t  =  t  +l/2At,  then  the  function 
F\(Uq,  . .  .LL  ,trj+l/2At)  used  in  the  Runge-Kutta  step  is  replaced  by 
Fj(S(Uii ,t+l/2/\t)  ,ui » •  •  •  5Uj,t+l/2At) .  Also,  at  the  end  of  the  step  the 
value  of  Uq  computed  by  the  ODE  solver  is  replaced  by  S(UQ,t)  provided 
Ug  is  still  in  the  constrained  set  Uj .  Obviously  this  requires  modifi¬ 
cation  of  the  ODE  solver.  There  is  no  guarantee  that  this  method  will 
converge.  In  fact,  as  we  will  see  shortly,  it  does  not  always  converge. 
The  algorithm  can  be  implemented  as  part  of  a  PDE  package  once  the  user 
has  supplied  the  subroutines  to  evaluate  p,  the  sets  Uj,  and  the  functions 
S^Ujj.t). 

The  second  method  is  a  generalization  of  the  differentiated 
boundary  conditions  described  in  section  3  in  equations  (9)  and  (10). 

In  this  case  the  user  is  allowed  to  reset  the  time  derivatives  used 
by  the  ODE  solver  to  compute  the  boundary  values,  that  is 


dU 


m,Q 


dt 


=  F„,,o%*£o*t)  and 


dU, 


m,J 


dt 


Fm,j  =  <16> 


Here  we  assume  a  system  of  equations  for  the  unknowns  u  •  where  there  are 

m,j 

M  unknowns  (1 <m<M)  and  the  mesh  points  are  X .(a=xn<xn <  ...<x1=b). 

j  u  i  u 
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T 

The  vectors  Fg  =  (F^  q,...,Fm  g)  and  F^  are  the  time  derivatives 
obtained  using  one  sided  difference  approximations  in  the  hyperbolic 
system  at  the  boundary.  The  vectors  Ug,Fg  and  the  time  t  are  supplied 
to  a  user  written  subroutine  which  must  then  determine  the  set  I  and 
return  values  for  - -jjp  =  F^  g,  for  pel.  The  remaining  time  derivatives 
for  p£l  are  the  values  F^  Q  obtained  from  onesided  differences  in  the 
hyperbolic  system.  This  method  is  going  to  be  difficult  to  explain  to  a 
user.  However,  it  is  the  only  method  that  has,  so  far,  worked  reliably. 

We  will  illustrate  this  method  by  the  following  example.  This 
is  a  system  with  characteristics  of  different  sign.  Chu  [3]  and  Sundstrom 
[10]  have  noted  the  difficulties  of  setting  boundary  conditions  for  such 
systems. 


Thi 


Ou  i 

at 

3  3u-j 
3x~ 

4  3u0 

-  »x 

0«t 

(17) 

au2 

29u-j 

33u2 

0<x<b 

3t 

"1)T 

3x 

system  is 

derived 

from 

3u  . 
at  ' 

_3u 

3x 

U  =  U-j  “Up 

u-j  =2u+v 

3v 

dy 

v  =  2u2“U-j 

U2=u+v 

at  : 

3x 

Therefore  the  following  boundary  conditions  are  proper,  since  they 
amount  to  a  specification  of  the  characteristic  variables  on  the 
inflow  boundary. 

at  x=0  v=2u2>(0,t)-u.j  (0,t)=-sin2Trt 

at  x=b  u=Ui  (b,t)-u2(b,t)=sin2ir(b+t) 


(18) 
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We  have  chosen  the  boundary  conditions  to  correspond  to  the  following 
solution 

u-,(x,t)  =  2sin27r(x9t)+sin27r(x-t) 

(19) 

u2(x,t)  =  sin2Tr.(x+t)+sin27r(x-t) 

To  use  the  first  method  of  setting  the  boundary  conditions  we  must 
specify  the  set  Uj  at  each  boundary  point.  There  is  no  unique  choice 
here,  since  neither  u-j  nor  u2  are  characteristic  variables.  We  will  try 
to  specify  u-|  at  each  boundary  from  the  given  boundary  conditions,  namely 

at  x=0  u-|  (0,t)=2u2(0,t)+sin27Tt  j20^ 

at  x=b  Ui (b,t)=u2(b,t)+sin2n(b+t) 


In  this  case  I  =  {!},  Uj  =  {u-j},  II  =  {2},  Ujj.  =  (u2>,  and  p  =  1  at 
both  boundary  points. 

A  derivative  rather  than  a  constrained  boundary  condition  can  be 
obtained  by  differentiation  of  the  above  equation,  namely 


+  2'rrCOs2'iTt 


+  2nCOS27r(b+t) 


at  x=0 

at  x=b 


(21) 


The  derivative  du2/dt  on  the  right  can  be  computed  from  the  hyperbolic 
system  using  one  sided  differences  and  then  used  in  the  user  supplied 
routine  to  compute  du-j/dt  by  equation  (21)  above. 

As  our  results  show  neither  of  these  methods  given  by  equations 
(20)  and  (21)  work  sati sfactori ly.  They  both  specify  the  inflow  character¬ 
istic  variable.  The  outflow  characteristic  should  be  computed  using 
one  sided  differences.  In  equation  (20)  the  inflow  characteristic  is 
specified  by  the  boundary  constraint.  However,  there  is  certainly  error 
in  the  computed  value  of  u2  used  on  the  right  side  of  the  boundary 
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constraint.  This  error  can  be  transmitted  to  the  other  boundary  and 
reflected  back.  The  boundary  condition  probably  should  not  allow  much 
error  in  the  incoming  characteristic. 

We  tried  a  third  type  of  boundary  condition  obtained  by  differen¬ 
tiating  the  boundary  constraint  and  combining  it  with  the  equation  for  the 
outgoing  characteristic  variable  obtained  from  the  hyperbolic  system. 

That  is,  at  x  =  0 

dul  0  du2  0 
- L’-  +  2-4^  =  -2'rrCOS2irt 


dt 

du 


1  ,0 


dt 

du 


2,0 


dt  dt 
Here  F-j  Q  is  an  approximation  to 
3-&U-J  4,3  u2 


F1,0"F2,0  * 


3X 


BX 


and  F 


2,0 


^x  $x 

obtained  using  one  sided  differences.  These  equations  yield 

du- 
dt 
du. 


'1,0 


2F1.0-2F2.0-2lrt:os2rt 


(22) 


2,0  _ 


dt 


F1  0_F2  q— 27rCOS2Trt 


There  are  errors  in  computing  F-j  g  and  F^  Q,  but  these  will  cancel  out 
in  the  computation  of  the  time  derivative  of  the  inflow  characteristic 
(v=2u2“U-j)  when  this  boundary  condition  is  used.  Perhaps  this  is  the 
reason  for  the  superior  performance  of  condition  (22)  over  (20)  and  (21), 
However,  we  do  not  have  a  solid  theoretical  understanding  of  these  results. 
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5.  Some  computational  results.  These  results  all  refer  to  the 
solution  of  equation  (17)  using  a  fourth  order  centered  finite  difference 
approximation  in  the  interior  and  third  order  one  sided  differences 
near  the  boundary  to  approximate  the  spatial  derivation  3/3 x.  The 
Runge-Kutta-Fehl berg  [5]  method  was  modified  to  allow  use  of  the 
"constrained"  boundary  condition  (20).  The  "derivative-constrained" 
condition  (21)  and  the  "derivative-characteristic"  condition  (22)  were 
also  used.  The  parameter  e  refers  to  the  error  tolerance  used  in  the 
Runge-Kutta-Fehl berg.  The  variable  J  is  the  number  of  mesh  points,  and 
x=b  is  the  right  boundary.  The  results  depend  on  b,  probably  because 
of  the  way  the  error  is  reflected  between  the  two  boundaries.  The 
error  is  the  relative  error  in  the  computed  solution  at  the  indicated 
time  t=T.  The  parameter  is  the  number  of  evaluations  of  the  time 
derivative  required  in  the  integration.  Each  time  step  requires  6 
evaluations  (5  if  it  follows  an  unsuccessful  step). 

There  seems  to  be  little  difference  between  the  results  for  the 
constrained-boundary  (20)  and  the  derivative-constrained  method  (21), 
except  for  a  slight  difference  in  the  number  of  functional  evaluations. 
This  difference  can  be  largely  eliminated  by  omitting  the  error  estimate 
for  the  constrained  boundary  variables  -  at  least  this  was  our  experience 
for  the  single  equation  (1).  Only  the  characteristic  derivative  method 
(22)  is  free  from  the  error  growth  which  is  probably  due  to  multiple 
reflections  from  the  boundaries.  Note  that  the  severity  of  the  error 
growth  depends  on  the  length  of  the  interval  (the  parameter  b).  The 
error  reinforcement  upon  reflection  is  probably  dependent  on  the  phase 
angle  which  in  turn  depends  on  b.  Of  course,  these  results  are  based  on 
a  single,  simple  test  case  and  may  not  apply  to  a  given  problem. 
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These  computations  were  performed  on  the  CDC  6400  at  the 
University  of  Colorado. 
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Table  I.  Behavior  of  the  matrix  A 

of  the  semidiscrete  scheme  u.'=Au_+£. 
Here  ^  denotes  the  real  part  of  an 
eigenvalue  of  A. 
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Table  II  .  Error  for  various  schemes 
applied  to  equation  (16) 
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Table  III.  Error  for  the  solution  of 

equation  (21).  Here  T=time, 
and  (lull  is  the  maximum  norm  of 
the  solution. 
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Table  IV.  Error  for  solution  of  equation  (17) 
with  various  boundary  conditions. 


